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With the on going construction of several large and medium scale laser in- 
Qf-^ • terferometric gravitational wave antennas around the globe, the problem of 
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<^ • the detection of gravitational waves has acquired great impetus. Since grav- 
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itational wave signals from astrophysical sources are expected to be weak, 
despite the state of the art technology being employed, the development of 
optimal signal extraction techniques and the consequent accurate determina- 
tion of the parameters of the signal is of major importance. Coalescing binary 



^ . systems are one of the most promising sources of gravitational waves. One 



reason is that such sources are easier to model and thus one can design de- 
tection strategies tuned to such signals. A lot of attention has been devoted 
in the literature studying such techniques and most of the work has revolved 
around matched filtering and maximum likelihood estimation. 

In a previous work, Monte Carlo simulations were carried out of the de- 
tection process using matched filtering for the initial LIGO/VIRGO config- 
uration for the first post-Newtonian corrected coalescing binary waveform. 
We had compared the results of our simulations with available estimates, 
obtained from covariance matrix considerations, of the errors in the deter- 
mination of the parameters. Our results showed that the covariance matrix 



underestimates, by over a factor of two, the actual errors in the estimation 
of parameters even when the signal-to-noise ratio (SNR) is as high as 10. 
Sources having SNR higher than 10 are expected to be few and hence this 
issue is of major concern. 

In this paper we probe the question as to why the Monte Carlo simula- 
tions give such high errors as opposed to those obtained via the covariance 
matrix. We present, a computationally viable statistical model of the dis- 
tribution, of the maximum likelihood estimates (MLE), of the parameters. 
This model reproduces the essential features of the Monte Carlo simulations, 
thereby explaining the large root mean square errors in the estimates, ob- 
tained in numerical experiments. The chief reason for the large errors seems 
to be the fact that the probability distribution of the estimated parameters 
is multimodal. Though only the central peak (corresponding to the actual 
values of the parameters) is dominant, the subsidary peaks occur 'far' away 
thereby contributing to large variances. We suggest that the variance or the 
standard deviation of an estimated parameter may not provide the best mea- 
sure of the error, for the kind of situation we encounter here. We therefore 
propose another criterion by which the MLE should be judged. 

In order to illustrate the model we have considered the Newtonian as well 

as the first post-Newtonian corrected waveform. We have assumed Gaussian 

noise, with a power spectrum typical of the LIGO/VIRGO type of detectors. 

The model we have used, however, is quite general, and robust, and will be 

relevant to many other parameter estimation problems. 
PACS numbers: 04.30.-Fx, OAM.+z 
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I. INTRODUCTION 



Large scale laser interferometric detectors of gravitational waves, namely, the LIGO 
IjlJ and VIRGO and the medium scale detectors, GEO and TAMA are expected to be 
operational by the turn of this century. Compact coalescing binary systems of blackholes 
and/or neutron stars are relatively 'clean' systems to model during their inspiral and their 
inspiral waveform can be predicted with a fair degree of reliability. This makes them the most 
promising sources for broad band detectors, and in particular, the upcoming interferometric 
detectors cited above. Binary systems are also valuable sources of astrophysical information 
as one can probe the universe up to cosmological distances. For instance, statistical analysis 
of several binary coalescences enables the estimation of the Hubble constant to an accuracy 
better than 10% Events that produce high signal-to-noise ratio can be potentially used 

to observe such non-linear effects, as gravitational wave tails, and to put general relativity 
into test in the strongly non-linear regime 0. Due to the weak coupling of gravitational 
radiation with matter the signal waveform has a very low amplitude and will not stand above 
the detector noise. In addition to the on-going efforts to reduce the noise, and hence increase 
the sensitivity of the detector, a considerable amount of research activity has gone into the 
development of efficient and robust data analysis techniques to extract signals buried in very 
noisy data. For a review on gravitational waves from compact objects and their detection 
see Thorne 

Various data analysis schemes have been proposed for the detection of the 'chirp' wave- 
form from such systems. Among them the technique of matched filtering is the most promis- 
ing ||9|-[TI|. Briefiy, this technique involves correlating the detector output with a set of tem- 
plates, each of which is tuned to detect the signal with a particular set of parameters. In 
order to obtain a high value for the correlation the signal waveform should be known to a 
high level of accuracy. The matched filtering technique is very sensitive to the phase of the 
signal and even if the template and the signal mismatch by even half a cycle the correlation 
integral is drastically reduced. 
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The fully general relativistic waveform from a coalescing binary system of stars is as yet 
unavailable. In the absence of such an exact solution, there have been efforts to find solutions 
perturbatively. Most of the work in this area strives towards computing the waveform correct 
to a high degree of accuracy so that the theoretical templates based on this will obtain the 
maximum possible SNR. 

The signal is said to be detected, if the maximum value of the correlation over all the 
parameters of the signal crosses a preassigned threshold which has been set by the false 
alarm one is prepared to tolerate. Once the signal is detected, the maximum likelihood 
estimates (MLE) of the parameters of the binary are those of the template with which the 
maximum correlation is obtained. The errors involved in such an estimation have been 
worked out by several authors P,p!^-[T8|, for the case of 'high' SNR and for the Newtonian 



and post-Newtonian waveforms using a single and a network of detectors. In |jT9[ exhaustive 
Monte Carlo numerical simulations were carried out to compute the errors in the estimation 
of parameters and covariances among them, for the case of the initial LIGO configuration 
taking into account only the first post-Newtonian corrections and assuming circular orbits. 
It was found that the errors as obtained from the simulations were larger by a factor of two 
or more from the errors as computed via the covariance matrix at astrophysically relevant 
SNRs. The discrepancy disappears as the SNR increases beyond certain value - typically 
15 in the Newtonian case and 25 in the post-Newtonian case. The comparision with other 
less stringent lower bounds has also not resolved the discrepancy. Nicholson and Vecchio 
pO| have recently computed Bayesian bounds such as the Weiss- Wainstein and Ziv-Zakai, 
for the case of Newtonian signals. They conclude that though these bounds are tighter than 
the Cramer-Rao bound, the numerical errors are still much larger than the computed lower 
bounds. 

In this paper we explain the discrepancy between the Monte Carlo simulations and 
the results obtained from the covariance matrix. We demonstrate that the probability 
distribution of the estimated parameters cannot be got from 'local' methods such as the 



ones used earlier . Since our main purpose in this paper is to demonstrate the validity of 



the statistical model, we initially use the simplified Newtonian model of the waveform. In the 
Newtonian case there are fewer parameters of the signal that one has to deal with, making 
the investigations simpler, analytically as well as computationally. We then specialize to the 
first post-Newtonian case. Here, although the discrepancy is larger, the problem is similar 
at a qualitative level. 



Following Finn [|12| we obtain an equation which relates the estimated parameters to the 
actual parameters of the signal for a given noise realisation. This equation is non-linear. 
Finn linearises the equation and obtains the errors in terms of the covariance matrix. Here 
we do not linearise the equation. We find that the equation has multiple solutions for 
the parameters in a certain sense and it is this multiplicity of roots forming islands in the 
parameter space which contributes significantly to the errors. Thus the problem is of a 
global nature and any local approximation will be inadequate in explaining away the errors. 
Moreover we suggest that the variance/covariance of the parameters is not a proper measure 
of the errors in this case. We therefore propose a new criterion by which the MLE should 
be judged. 

The paper is organized as follows. In section ^ we briefiy describe the gravitational 
wave signal waveform and the parameters on which it depends, namely, the amplitude, the 
time of arrival, the phase at arrival and the 'chirp times', which characterizes the time 
for the binary to evolve from some fiducial time to the actual merger. These parameters 
are found to be very convenient when we carry out Monte Carlo simulations. It turns out 
that the covariance matrix is independent of these parameters and hence it is sufficient to 
carry out the simulations only for a particular set of parameters. Further in this section we 
also describe the characterstics of the noise that we assume will be present in the detector 
and briefiy review the matched filtering technique. In section |ITT| we present the results 
of the Monte Carlo simulations for the Newtonian case. We carry out transformations of 
the parameter space which bring out the chief features of the distribution of the estimated 
parameters. We show that the estimated parameters do not lie in a simply connected 
region around the signal parameters, but instead are distributed in multiple 'islands' in the 



parameter space considered. We next present a geometric representation of our statistical 
model. We then apply this model to the Newtonian chirp waveform, and compare the model 
with the Monte Carlo simulations. In section |rV| we deal with the post-Newtonian waveform. 
Finally in section |V| we summarise our results. We propose an alternative measure for the 
error which performs reasonably better than the variance as a measure of the error. 

II. THE SIGNAL AND THE NOISE 
A. The Chirp Signal 

When constructing templates for on-line detection, it is sufficient to work with the so 
called restricted post-Newtonian gravitational waveform. In this approximation the post- 
Newtonian corrections are incorporated only in the phase of the waveform, while ignor- 
ing corresponding corrections to the amplitude ||2l]|. Consequently, the restricted post- 
Newtonian waveforms only contain the dominant frequency equal to twice the orbital fre- 
quency of the binary computed up to the relevant order. In the restricted post-Newtonian 
approximation the gravitational waves from a binary system of stars, modeled as point 
masses orbiting about each other in a circular orbit, induce a strain s{t) at the detector 
given by 

s{t) =A{nf{t)y/' cos [^{t) + ^, (2.1) 

where f{t) is the instantaneous gravitational wave frequency, the constant A involves the 
distance to the binary, its reduced and total mass, and the antenna pattern of the detector 
and $ is the initial phase of the wave at some fiducial time t = tg. The phase of the waveform 
ip{t) contains several pieces corresponding to different post-Newtonian contributions which 
can be schematically written as 

^{t) = ipoit) + ifiit) + ypi.5(t) + . . . . (2.2) 

The evolution of the phase depends on the masses of the two components of the binary, 
characterised by the reduced mass, fi and the total mass, M of the system. Here <fo{t) is the 
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dominant Newtonian part of the phase and ipn{t) represents the nth order post-Newtonian 
correction to it. The Newtonian part of the phase is sensitive only to a particular combination 
of the masses of the two stars, frequently characterised by its 'chirp mass', 
We give below the waveform correct to the first post-Newtonian order: 



Mi) 



167r/,ro 



^pl{t) = 4:71 fsTi 

where f(t) is given implicitly by. 
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where tq and ri are constants having dimensions of time given by 
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where rj = /i/M, and fs is the instantaneous gravitational wave frequency of the signal at 
t = tg. The time elapsed starting from an epoch when the gravitational wave frequency is 
fs till the epoch when it becomes infinite will be referred to as the chirp time of the signal. 
In the quadrupole approximation tq is the chirp time whereas it is tq -|- ti for the first post- 
Newtonian case . The Newtonian part of the phase is characterised by three parameters: 
(i) the time of arrival tg when the signal first becomes visible in the detector, (ii) the phase $ 
of the signal at the time of arrival and (iii) the chirp mass. At this level of approximation two 
coalescing binary signals of the same chirp mass but of different sets of individual masses 
would be degenerate and thus exhibit exactly the same time evolution. This degeneracy 
is removed when post-Newtonian corrections are included. The parameters tg and $ are 
kinematical that fix the origin of the measurement of time and phase, respectively, while the 
Newtonian and the post-Newtonian chirp times are dynamical parameters in the sense that 
they dictate the evolution of the phase and the amplitude of the signal. 



The parameters tq, ti and ts have the dimensions of time. We convert them to dimen- 
sionless parameters by multiplying them with 27!- fg. Thus we have the parameter set, 



At = n^] = {A, 2TTfsts, $, 27r/,ro, 27r/,ri}. 

In the stationary phase approximation the Fourier transform of the restricted second- 
post-Newtonian chirp waveform for positive frequencies is given by . 



s{f) = yW/-^/^exp 



7r 



(2.6) 



where A, is the amplitude parameter depending on the distance to the binary, as well as the 
chirp mass and A/" is a normalization constant to be fixed by means of a scalar product to 
be introduced later, and 



Xi 



I 

X2 = -1, 

fs 5^5 [fs 

X4 = y-2 + ^. 

Is J 



-5/3 



(2.7) 



For / < the Fourier transform is computed using the identity s{—f) = s*{f) obeyed by 
real functions s(t). 



B. The noise 



The output of a gravitational wave detector such as the LIGO, will comprise of data 
segments, each of duration T seconds, uniformly sampled with a sampling interval of A, 
giving the number of samples in a single data train to be = T/A. Each data train can 
be considered as a A^-tuple x = {x^, x^, . . . , x'^~^}, being the value of the output of the 
detector at time aA. The set of all such A^-tuples constitutes an A^-dimensional vector space 
V where the addition of two vectors is accomplished by the addition of corresponding time 
samples. For later convenience we allow each sample to take complex values. A natural basis 

8 



for this vector space is the time basis e° = where r and a are the vector and component 
indices respectively. Another basis which we shall use extensively is the Fourier basis. 

A gravitational wave signal from a coalescing binary system can be characterised by a 
set of parameters /x", a = 0, 1, ...,m — 1 belonging to some open set of the m-dimensional 
real space R"^. The set of such signals s{t; ^i) constitutes a m-dimensional manifold S which 
is embedded in the vector space V. Note that Greek characters in boldface denote the full 
collection of parameters characterising the signal. The parameters of the binary can be 
regarded as coordinates on the manifold. The basic problem of signal analysis is thus to 
determine whether the detector output vector x is the sum of a signal vector and a noise 
vector, X = s + n, or just the noise vector, x = n, and furthermore to identify which 
particular signal vector, among all possible. The latter is relevant to this paper where we 
are interested in estimating the parameters and also the errors made in such a estimation. 
Errors in the estimation arise because the noise contaminates the data. 

The noise in ground based laser interferometric detectors will have, in general, both a 
Gaussian and a non-Gaussian component. The main sources for the Gaussian component 
are the shot noise due to photon counting, the thermal noise in the mirror suspensions along- 
with the mirror itself and seismic noise. The non-Gaussian component can be contributed 
by numerous sources like sudden strain releases in the mirror suspension wires or even if 
lightning strikes. It should be possible to remove most of the non-Gaussian component by 
using environmental monitors and looking for coincidence among detectors located at widely 
separated sites. It is, therefore, assumed usually that the detector noise will be a Gaussian 
random process. Over a time scale of hours, it can also be assumed to be stationary. 

The power spectral density of the Gaussian noise component rises very steeply towards 
the low frequency end due to seismic effects. At the high frequency end it is dominated by 
photon shot noise which leads to a rise towards higher frequencies. Thus the data will have 
to be bandpassed with a low frequency seismic cutoff, /s, and a high frequency cutoff, fc- 
We use the power spectral density expected for the initial LIGO as given in [||]. Accordingly, 
we choose fs = 40 Hz and fc = 800 Hz. 
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In the absence of the signal the output will contain only noise drawn from a stochastic 
process which can be described by a probability distribution on the vector space V. We 
assume that the noise has its mean zero, or that is, ri" = 0, where the overbar denotes an 
ensemble average. Then the covariance matrix of the noise C"'' is defined as, 



C 



ab 



arnb 



n"-n 



(2.8) 



If the noise is assumed to be stationary and ergodic then there exists a noise autocorrelation 
function K{t) such that C"* = K{\a — b\A). In the Fourier basis it can be shown that the 
components of the noise vector are statistically independent |]TU[ and the covariance matrix 
in the Fourier basis will contain only diagonal terms whose values will be strictly positive: 
Qaa _ fiafj^*a_ xhis implies that the covariance matrix has strictly positive eigenvalues. The 
diagonal elements of this matrix C"*^ constitute the discrete representation of the power 
spectrum of the noise 

Gaussian noise can be described by the distribution. 



Po{n) = MnBxp 
where M„ is a normalization constant given by 



N-1 



7^ E [C-%bn^n' 

a,b=0 



(2.9) 
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Equivalently in the Fourier domain this can be written as, 

1 N-1 



Po(n) = M„exp 
= M„exp 



n 



a,b=0 
N-1 



a=0 



(2.10) 



where in the last step we have used the diagonal property of the matrix C"^ which implies 
that [C-%a = 1/C°". 

In the presence of the signal s(/i) the above probability density function (pdf) gets 
modified but in a very simple way since we have assumed that the noise is additive. We 
have, 
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Pi(x) = po(x - s(/i)), 



(2.11) 



where pi(x) is the pdf of x when the signal s(/i) is present in the data stream. 

In this paper we shall assume the noise to have a power spectrum consistent with the 



initial LIGO instrument. We use the fit to the noise curve as given in ||16 

-4 / 
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(2.12) 



The value of 5*0 will not concern us since what matters is only the ratio of the signal amplitude 
to that of the noise, in other words the SNR. We shall set the value of Sq to be unity and 
accordingly adjust the amplitude of the signal to get the required SNR. 



C. The Matched Filter 

In the absence of any prior information it must be assumed that all the parameter val- 
ues, within their respective range, are equally likely to occur. In such a case, the method of 
maximum liklihood can be used. When the noise is a stationary Gaussian random process, 
the method of MLE reduces to the so called matched filtering technique. Matched filtering 
involves correlating the detector output with a bank of matched filters each of which corre- 
sponds to the signal waveform for a fixed set of parameters. To this end, we define a scalar 
product on V. In the continuum limit the scalar product between two vectors x and y is 
given by, 

(x, y) = r df ( x{f) r (/) + x*{f) y{f) ) , (2.13) 

where, the Hermitian property of the Fourier transform of a real function has been used 
to restrict the domain of integration to positive frequencies. Sn{f) is the power spectral 
density of the noise. The Fourier domain is convenient since stationarity of the noise has 
been assumed. The norm of a vector z will be denoted by ||z|| = (z,z) . In eqn. ( |2.6| ) we 
had left J\f undefined. We choose the value of A/" such that ||s(/x)|| = A. From the definition 
of the scalar product in eqn. ( |2.13|) and from eqn. (|2.6| ) it follows that. 
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(2.14) 



The integrand in the above equation, /(/) is plotted in Fig. 0. This function peaks at around 
/ = 135Hz and a major contribution to the integral comes from a region around 135Hz. 

We shall use normalized matched filters, that is,we set h(/i-') = s{ii)/A. Henceforth we 
shall use the symbols a,b, . . . to indicate indices whose range of values includes which as 
a index denotes the amplitude parameter, i.e. fi^ = A. Indices . . . will not include the 
amplitude parameter and will never take the value 0. 

The output called the correlation c(/i-') of the matched filter with parameters fi^ is then 
just the scalar product given by. 



Given the data vector x, the correlation is computed for the entire feasible range of param- 
eters, continuously over the kinematic parameters tg and $ and for discrete values of the 
dynamical parameters tq and ti. The filters corresponding to the discrete values of tq and 
Ti constitute the filter bank. The maximum of c{fi^) is computed and compared with a pre- 
assigned threshold determined from the false alarm probability. A detection is announced 
if the maximum of c{fi^) crosses the threshold. The parameters fi^ for which the correlation 
is maximised are the MLE of the parameters of the signal. However, these will in general 
differ from the actual signal parameters fi^ due to the presence of noise. The difference 
A/i-' = fi^ — fi^ is the error made in estimating the parameters. When the signal is weak as 
compared to the noise (low SNR), the estimated parameters in general will differ by a large 
margin from the true ones, while in the limit of high SNR, the two will almost coincide. 
Thus in general Afi^ is not small. 
Let the data vector be 




(2.15) 



x = Ah.{fl^) + n. 



(2.16) 



where A is the amplitude of the signal and fi^ are the remaining parameters. We have 
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separated the amplitude as it occurs in the signal as just a scale factor multiplying the 
signal. The SNR is defined as, 

p={s{il),sifi)f' = A. (2.17) 



From eqs. (|2.6| ) and (|2.13|) we see that the SNR is equal to the amplitude parameter of the 
signal, A. 

When we consider an ensemble of noise realisations, the c{fi^) becomes a random variable. 
For a fixed set of parameters /i-^, c(/i-') is a Gaussian random variable since it is a linear 
function of the noise (eqn. (|2.15|) . The process of maximization over the parameters is 



however a nonlinear process and hence the statistics of c{fi^) will be non-Gaussian. For a 
given realization of noise, we assume that the global maximum of c{fi^) will also be a local 
maximum. This assumption depends on the specific nature of the waveform and the filter 
bank. The equations which will be satisfied at the estimated parameters fi^ = fi^ , are, 

dc 

In the limit of high SNR the mean square errors in the measurement of the parameters 
are characterised by the so called covariance matrix C"*'' [|10| defined as. 



c 



ab 



ab 



(2.20) 



From the expression for the Fourier transform in the stationary phase approximation 
(eqn. (p.6|)) and from the definition of the scalar product (eqn. ( |2.13| )) we see that the phase 
function in the Fourier transform simply cancel out and only the amplitude remains. We 
also find that oc A'"^. 

The idea now is to compare the results of the Monte Carlo simulations with those ob- 
tained via eqn. ( p.l9| ). In the limit of high SNR, it is expected that the errors will agree with 
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those obtained from the covariance matrix. In the following sections we briefly mention the 
results of the Monte Carlo simulations and then analyse in detail eqn. ( p.l8| ). The goal is to 
first check the agreement between the two and then understand how the additional errors 
arise by studying the consequences of eqn. ( p.l8| ). 



III. THE NEWTONIAN WAVEFORM 



A. Monte Carlo Simulations 



In this section we shall restrict ourselves to the Newtonian waveform. We have the four 
parameters: 

= /i" = ^\ /x^} = {A, 27r/A, 27r/,ro}. 

Simulations were carried out for a test case of tq = 5.558 s This corresponds to a binary 
comprised of a lOM© black hole and a IAMq neutron star. The waveform was cutoff at a 
frequency of 800 Hz and the data train was sampled at 2000IIz. As mentioned before the 
power spectral density was chosen to be consistent with the initial LIGO interferometer. 
The range of tq for the filters was [5.358, 5.758] with a spacing of 1 msec. The simulations 
presented here were carried out for an SNR of 10, i.e., A = 10. We considered 12000 
realizations of noise. 

The covariance matrix for the Newtonian case is given below: 



C 



(m) 



(m) 



1 

l2 



0.0 0.0 0.0 

0.0 222.50 322.26 -227.25 
0.0 322.26 469.16 -328.65 
0.0 -227.25 -328.65 232.26 



(3.1) 



In the above the diagonal values denote the variances in the parameters and the nondiagonal 
components denote the covariances between the parameters. Since the components do not 



depend on the parameters , (see p2| , |T9| for details) we can choose a new set of parameters 
in which the covariance matrix is diagonal. The transformation can be affected by means of 
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an orthogonal transformation. We have the new parameter set v = {A, u^, u"^, z/^}, related to 
fj, by means of the relation u = Sn, where S is a orthogonal matrix, which for our particular 
case is: 

'^1.0 0.0 0.0 0.0 ^ 



S 



0.0 -0.794 0.129 -0.594 
0.0 0.359 -0.690 -0.629 
0.0 -0.491 -0.713 0.501 



(3.2) 



In this new coordinate system the covariance matrix C(i,) = SC(^i)S ^, and is given by, 



C(^) 



i2 



^ 0.0 0.0 0.0 ^ 



0.0 0.024 0.0 0.0 
0.0 0.0 1.640 0.0 
0.0 0.0 0.0 922.3 



(3.3) 



In the high SNR limit the root mean square errors in the parameters are given by. 



(J,, a 



(3.4) 



For an SNR of 10 the values are {1.0, 0.015, 0.128, 3.037}. 



In a previous work [|l^ we had performed detailed Monte Carlo simulations to study the 
variation of the errors with the SNR. The simulations carried out in that paper were with 
a slightly different power spectrum, as we have used a simple fit to the noise curve in this 
paper. We reproduce in Fig. ^ the variation of errors in the parameters with the SNR as 



given in Figure 5 of for the parameters tq and t^, except that we use the log scale for 
both the axis, as is conventional in statistical literature. The continuous line represents the 
errors computed via the covariance matrix and in this approximation the errors are inversely 
propotional to the SNR. The dotted line represents the errors as obtained from the Monte 
Carlo simulations. The rest of this paper tries to explain the discrepancy between the two 
curves in this figure at low SNRs p ~ 10. 
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B. The equation for the errors 



The expression for the Fourier transform of the chirp in eqn. ( |2.6| ) in the new coordinates 
retains its form, but the functions Xi{f) are now transformed to, 



3 



Rewriting eqn. ( |2.18| ) in the new parameters we have at = 0^ 




dh 

97 



(3.5) 



Using the definition of the scalar product in eqn. (|2.13|) , we have. 




(3.6) 



where Au^ = 0^ — In the Newtonian case these constitute a set of three nonlinear 
equations connecting z/* to Kj. The quatities are random variables. Since z>* depend upon 
n, are in general nonlinear functions of the noise. This implies that the statistics of /tj can 
in general be non-Gaussian. However, for chirp signals of the type we consider, the statistics 
of Ki and K2 will turn out to be Gaussian. This will be demonstrated in the next section. 
We define the total phase difference 9 between the signal and the filter as. 



which is relevant for future considerations. In the high SNR limit the errors Az/-' are small 
and hence we can make the approximation sin(0) ~ 9. This corresponds to using the 
covariance matrix to provide an estimate of the errors. At astrophysically interesting SNRs 
10 this assumption does not hold. 

It is clear that a good correlation between two waveforms is obtained when the phase 
difference between the two waveforms is a multiple of 2ti and is roughly constant over the 
time for which the waveforms last. This assumption is found to be true in the simulations. 
Even though there are far more outlying points than what is predicted by the covariance 



3 




(3.7) 
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matrix the phase difference between the signal and the optimally matching template is found 
to be roughly constant. A typical case is illustrated in Fig. ^, in which 9 is plotted as a 
function of frequency for Az/* = {1.62,-8.49,7.3}. The function 9{f) is found to be close 
to — 47r in the region around 135 Hz, from which the maximum contribution to the SNR is 
obtained. Note that the values of Az/^ and Az/^ are way beyond the rms errors predicted by 
the covariance matrix. 

Since the function /(/) in Fig. |^ peaks at / ~ 135Hz., we define Om to be the value of 9 
at / = 135Hz., i.e., 

9m = -3.87615Az/^ + 0.739347Ai/2 + -0.0149334Az/3. (3.8) 

In Fig. ^ we illustrate how the 9m is distributed across many realizations of noise. It is clear 
from the figure that 9m is strongly clustered around multiples of 2tt. Since 9m is a linear 
combination of the parameters, it is clear that at least some of the parameters will show 
similar clustering properties. The variable 9m is a very convenient indicator of the rough 
location of the parameters in the parameter space. If the absolute value of 9m is large then 
the estimated parameters are far away from the actual value of the parameters. 

Fig. ^ is a z/^ — z/^ scatter plot, where each point corresponds to a realization of noise. 
The plot shows that z/^ and z/^ parameters are clustered in well separated 'islands' in the 
parameter space considered. The variable 9m computed for points in any specific island 
yield values close to a specific integral multiple of 27i depending upon the island we choose. 
Multiple solutions to eqn. ( p.6|) are responsible for the islands. Since the pdf in the Kj space 
is largest at the origin, the islands correspond to 9{f) ~ 2kTT, where k is an integer. We 
give an explanation of this below. The phase parameter /i^ = $ is constrained to lie in the 
interval [— 7r,7r]. Using the relation /i* = [S^^^jU^ we have, 

$ = /i^ = 0.129212z/^ + -0.689524Z/2 + -0.712644z/l (3.9) 

Since the rms error in as calculated from the covariance matrix at an SNR of 10 is 3.04, 
it is highly likely that the absolute value of the last term in eqn. ( p.9|) , ( which dominates 

17 



the other terms in the same equation, since the rms errors for the other parameters are 
much lesser) exceeds vr. This forces the parameters z/^ and to jump to values such that 
$ remains in the range [— 7r,7r]. In order to calculate the amount by which the parameters 
jump, we consider eqn. (|37^). Since Qm must be a multiple of In in order to get a good 
match, each of the first two terms on the right hand side of eqn. ( |3.8| ) contributes vr. This 
implies that jumps by vr/S. 87615 ~ .81 and jumps by tt/O. 739347 ~ 4.25. Since the 
two terms are of opposite signs, we find that increases when decreases and vice versa. 

We now revert back to eqn. ( p.6|) . Since Q is close to an integral multiple of 27r in the 
frequency region of interest, we can write, 

sin(^^) ^ I ^r7,(/)Az/^- -2A;7r I , (3.10) 



for some integer k. We have. 



^^^^^ r f!MmnM,f _ r H^a. (s.n) 



Sn{f) Jo Snif) 



Since the u coordinate system is orthogonal. 



Snif) 



df = for t^j, (3.12) 



and 



Therefore, 

= 2A;7rG„ (3.14) 

where, 

G, = 2Af' rQ^df. (3.15) 

The amount by which the parameters z/^ and u"^ jump can be calculated more accurately 
using eqn. (|3.14|) . For two successive values of the integer k, the parameter has to jump 



by, 2tiGi/CI^ in order to have the same value of ni. This works out to be 0.812 for the 
parameter . Similarly the value for the parameter z/^ turns out to be 4.332. These values 
are consistent with the scatter plot in Fig. ^ 

It is clear from Fig. ^ that the distribution of the parameters and z/^ will be markedly 
multimodal. However, the distribution of is relatively smoother. This is illustrated in 
Fig. ^ which is a scatter plot on the v'^ — plane. Though there are gaps along the v'^ axis, 

takes all values in the range shown. The parameter 1/3 is relatively well behaved, i.e. it 
does not exhibit any sudden jumps like the other parameters. 

The variances as obtained from the Monte Carlo simulations in the parameters are, 
Sj,i = {0.421,2.26,2.53} whereas, the values predicted by the covariance matrix are a^i = 
{0.0153,0.1281,3.037}. In the case of the parameters and v'^ the observed variances are 
much larger than the lower bounds, due to the jumps which these parameters make. However 
we notice that the observed variance in the case of is actually less than the lower bound. 
This is due to the fact that the Cramer- Rao bounds are applicable only to parameters which 
are allowed to vary freely. For instance the variance for the parameter $ = /i^ can have 
the maximum value of vr^ whatever be the SNR. The restriction of the range of /i^ puts a 
constraint on the values of the parameters and this accounts for the observed error being 
less than the Cramer-Rao bound. In the next section we give a more quantitative model for 
the distribution of the parameters. 



C. Geometrical Perspective 

In this section we use differential geometry to arrive at a statistical model for the distri- 
bution of the parameters. The model described here is quite general and is applicable to the 
estimation of parameters obtained by means of the maximum likelihood method, for an arbi- 
trary signal in the presence of Gaussian noise. The set of signal vectors, s(i/) = s(t; i^), where 
v = {uq, . . . jUm-i}, is a m-dimensional parameter vector, will describe a m-dimensional 
manifold S embedded in V, the vector space of all detector outputs. (See [IT^P^ 
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introduction to the use of differential geometry in gravitational wave data analysis, and ||25 
for the application of differential geometry to statistics.) Let the output of a detector x, 
contain a signal s(i>). Then x = s(i>) + n, where n is a noise vector drawn from the noise 
distribution. The distance, D{iy) between x and a point s{i/) is given by, 

D{u) = ||x-s(i/)|| = (x-s(zy),x-s(l/))^/^ (3.16) 
= [(x, x) - 2 (x, sH) + (sH, sH)]^/2 . (3.17) 

The MLE of the parameters is that point £> on the parameter space which minimises 
D{y). This is equivalent to maximising the correlation c{i>) = (x, s(i^)) provided we keep 
{s{u),s{h')) constant. 

In the limit of high SNR, s(i>) will lie in a small region around s(i>) on the manifold, 
effectively the tangent space to the manifold at s(i>). In this case, the difference, s(/>) — s(i/) 
can be satisfactorily approximated as the projection of the noise vector onto the tangent 
space. This implies that s{l>) — s{i>) is linear function of n. Further, if the parameters form 
a Cartesian system of coordinates, then, they too will be linear in n and the distribution 



of the parameters can be described by a multivariate Gaussian [|T2[ . The covariance matrix 
of this distribution defines a lower bound on the errors in estimation and is termed as the 
Cramer-Rao bound. 

If the global minimum of D{u) is also a local minimum then, at u = u, 

dD{u)/du'' = 0, for a = 0,...,m-l (3.18) 

which implies, 

(s{i>)+n-s{O),^{P)'j = 0. (3.19) 

These are a set of m equations, one for each parameter. We interpret these equations 
geometrically as follows. The equations imply that the vector x — s(i>) is orthogonal to each 
of the coordinate basis vectors, d/du"", evaluated at u = u. Thus i> is a local extremum when 
the tip of the vector x lies on that N — m dimensional hyperplane Bo-, which passes through 
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s{u), and is orthogonal to the m-dimensional tangent plane at s(r>). This hyperplane, Bo, is 
the intersection of the (A^ — l)-dimensional hyperplanes, TV"?, each orthogonal to a coordinate 
basis vector d/dv°' at u = u. 

Let la be the normalized coordinate basis vectors at i>, i.e. 



1. 



ds 



ds 



(3.20) 



We define to be the minimal distance from s(i>) to the hyperplane A/"^, which is given by 



Ta = (S(Z>) - S{i>),\a) . 



(3.21) 



A schematic illustration of the above discussion is given in Fig. ^ 

The pdf for the vector x to lie on A/"^ can depend only on la and r^. If the vector x is to 
lie on Bp, then it must simultaneously lie on each of the normal hyperplanes, A/"^. The pdf 
for X to lie on Bo is given by the expression, 



Pira) 



'm—l 



n5(((n-rA),U) 



a=0 



p(n) d'^n, 



(3.22) 



where the S denotes the Dirac Delta function. Substituting for the Gaussian distribution for 
the noise p{n) in the equation above and integrating, we get, 



exp 



p(ro,r2, . . . ,r^_i) 



m—l 

-I E [7- 

a,b=0 



l-\ab 



[(27r)™det M^'/' 



(3.23) 



where, 7aft = (la,lfe)- The integration though tedious, is quite straightforward. Note that 
each of the Gaussian random variables will have unit variance as is obvious from the 
definition of 7ab. Moreover, the matrix 7^^ is very closely related to the Fisher information 
matrix Tab as defined in eqn. ( p.l9|) . Whereas the components of the Fisher information are 
got by taking scalar products of the tangent vectors on the manifold, the components of the 
matrix 7^6 are obtained by by taking scalar products of the normalized tangent vectors on 
the manifold. 
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D. Statistical model for the Newtonian Chirp 



We now specialise to the case of the Newtonian chirp. We use the parameters defined 
earlier, in the previous section. Since s(i>) = Ah{C'^), the above equations for la and give, 



In = Uu' 



h 



dh 



ro(z>^i)=i(h(z>^),h(^)'=))-i, 



dh 



(3.24) 
(3.25) 
(3.26) 
(3.27) 



Since the u coordinate system is orthogonal, 'jab turns out to be nothing but the identity 
matrix. The distribution of the variables is thus a joint Gaussian given by the expression. 



1 



27r' 



exp 



1 



a=0 



(3.28) 



The Tj variables are closely related to the variables defined in section [III B| in eqn. (p. 51): 



A 



dh 

du 



/TV' 

J- 77 



where we have used the definition of Tab in eqn. ( |2.19| ). Thus they differ only by a factor 
which is simply a constant from eqs. (p.l3|) and (p.6| ). This is a consequence of the intrinsic 



flatness of the manifold [jT9[ and the particular parameterization adopted. 



From eqn. (|3.28|) we would expect the marginal probability distribution of each of the 
variables to be a Gaussian distribution with unit variance. Using the ensemble of es- 
timated parameters from the Monte Carlo simulations we can obtain the histograms and 
consequently the distributions of the variables r^. In Fig. ^ we plot the probability distri- 
butions of the variables r^, and compare them with the expected Gaussian distributions. It 
is clear that though ri and r2 are Gaussian random variables to a good approximation, 
shows a pronounced non-Gaussian behaviour. The reason for this discrepancy can be traced 
to the fact that the phase parameter is constrained to lie in the range [— vr, tt]. In Fig. ^ we 
observe that in the central island the z/^ parameter gets abruptly cutoff at a value of about 
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4.5 and —4.5. Since the points on the central island are 'close' to the point corresponding 
to the actual value of the parameters, we can apply the eqn. ( |3.14| ) with k = 0. This gives 



a value of k,^ ~ 0.00488 and consequently ^ 1.5. We observe in Fig. ^ that the dip in 
the distribution of occurs at the same point i.e. at ^ 1.5. To further elaborate on 
this point, we plot in the Fig. |^, the three variables v/s the variable 6m defined earlier. 
The figure clearly illustrates that while ri and r2 take on there entire range of values in the 
central island, does not. This establishes a connection between the dip in the marginal 
distribution and the phase parameter being constrained to the range [— 7r,7r]. Although a 
deeper understanding of this is in order, we continue our analysis assuming the distribution 
of r3 to be given by a Gaussian. 

Had the map between i> to r been bijective, it would have been possible to write the 
distribution for the estimated parameters simply as a product of the pdf for the variables 
Ta and a Jacobian factor. However, a given set of values for r = {va} would in general 
correspond to multiple parameter vectors where the range of values of / depends on the 
number of solutions. This is clear from Fig. ^ where we observe that the same value of re- 
occurs for different values of the variable 6m, and hence in different 'islands'. We also observe 
that r3 takes almost all its possible values in the central island and the two adjoining ones. 
Moreover, we notice that it is approximately true that takes different values in each of 
the three islands. Thus if we restrict ourselves only to the central and two adjoining islands, 
then the map between and rj is bijective to a good degree of approximation. For a fixed 
set of values for {rj} we have therefore, a unique solution satisfying eqs. (|3.27| ), if we 



restrict ourselves to the three islands identified above. We shall henceforth term the islands 
identified above as the primary group of islands and the solution there will be called the 
primary solution. (It is to be emphasized that the above discussion is applicable only to the 
Newtonian waveform. As we shall see later for the post-Newtonian case, the primary group 
of islands contains more islands on either side of the central island.) There will of course 
be other solutions in the other islands. It is to be noted that number of multiple solutions 
for a given set of values for {rj}, depends on {rj}, and we do not imply that each island 
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admits a solution for a fixed set of values for {vj}. The reason for the term primary 
solution is to be found in Fig. |1^, where we have plotted a histogram of the variable 6m 
which gives us an idea as to how many points occur in each island. We observe that 99% 
of the points lie in the central and the two adjoining islands for the Newtonian case. Thus 
there is an overwhelming probability for the points to lie in one of these islands, and the 
primary solutions occur much more frequently as compared to the other solutions for a fixed 
value of rj. 

The problem now is to determine the various solutions for a fixed set of values of 
{tq}, and compute the probability that a particular u'^'^^ will be selected amongst others 
for a given r. This is essentially the probability that the amplitude A^™-) is greater than the 
amplitude at the other solutions. We shall denote this probability as P{i>'^'^^\v). For a fixed 
set of values of {r^}, we will have multiple solutions, {C'^Y^\ . . . to eqn. ( |3.27| ). The 

correlation obtained at these points will be 

i(') = (x,h(i>^-)«). 

It is to be noted that for a fixed {z^-* }, A will be a Gaussian random variable with a variance of 
unity. In order to calculate P{{i'^}^^''\{rj}) we need to identify all the solutions corresponding 
to a fixed set of values r^. The identification of the multiple roots is quite a problem, and so 
we make the following approximation. We assume that for one of the solutions, that is the 
primary solution which we shall denote by corresponding to {^j}, the probability 

P({z>-'}'^^^|{rj}) is almost unity. If this is true, then, we only need to compute the probability 
that the correlation at an arbitrary point on the parameter space exceeds the correlation at 
the primary solution point which shares the same set of values of {rj}. This corresponds to 
computing the probability of the Gaussian random variable A^^^ , exceeding A^^^ . Of course, 
P({i>-'}''^^|{rj}) is set to unity. The justification for the above procedure is essentially the 
fact that nearly 99% of the points lie in the primary group of islands. 

So for evaluating the distribution of the estimated parameters at = we follow the 
following procedure: 
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1. Determine rj{i'^) using eqn. (|3.27|) . 



2. Determine {z/*^}*^^^ using eqn. ( |3.14[ ) for an appropriate value of A;, i.e. k takes one of 
the values {—1, 0, 1}. 



3. Determine the probablility for A = {^x, h(z/'^)y to be greater than A*^') = <^x, h(z/'^)'^^^y. 
(If z>^ is already a primary solution then this probability is set to unity.) 

4. Set P{{i^"'}\{ra}) equal to the calculated probability. 

5. Write the the pdf of the estimated parameters as 



p(i>) 



(2vr)- 



; J {€>) P (i^|r) exp 



^ a=0 



(3.29) 



where, Jii^) is the Jacobian of the transformation from to the variables r^, which 
is essentially, 



J(C) = det 



(3.30) 



Since the amplitude parameter is not of primary interest to us, we shall integrate the 
distribution over A. We use a threshhold of 7, which means that we reject any realization 
whose measured value A is less than 7 in the simulations. The amplitude parameter enters 
only via tq. So the distribution for the remaining three parameters 0^ is. 



u |r I exp 



E 

i=l 



X 



(2vr) 



1/2 



7.0 



exp 



dA. 



(3.31) 



Since the SNR is chosen to be 10 the second factor in the equation above is very close to 
unity. 



We now compare the Monte Carlo results with the distribution given in eqn. ( 3.31| ). We 
compare the one dimensional marginal distribution p(z>^) with the histogram obtained via the 
Monte Carlo method. The marginal distribution p(z>^) is obtained by integrating eqn. (|3.31|) 
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with respect to and . This is done numerically. Though the parameter has the least 
root mean square error of .015 as predicted by the covariance matrix, its distribution has 



the most pronounced non-Gaussian behaviour. In plot (a) and (b) of Fig. |Tl] we display the 
distributions obtained from the Monte Carlo simulations and the statistical model 

respectively. Plots (c) and (d) in the same figure zoom in on the first two maxima occurring 
on the right of the central maximum. It can be seen that in the case of plot (d) the match is 
not very good even though the location of the peaks match fairly well. The difference here 
could come from either the Monte Carlo method or the statistical model. 



In the histogram of the variable Qm illustrated in Fig. |TO|we have seen that about 79.5% 
of the points in the simulations fall in the central island and about 10% in each of the 
adjoining islands at an SNR of 10. We can obtain the corresponding numbers from our 
theoretical model as given in eqn. ( |3.29|) , by integrating the distribution over all the param- 
eters in the domain corresponding to each island. We obtain the values 82% for the central 
island and about 9.5% for each of the adjoining islands. Thus the statistical model shows 
good agreement with the Monte Carlo results by this criterion. The number of points in 
the subsequent islands falls off rapidly. The contribution to the variance from each island 
depends on the number of points in that region and the location of island in the parameter 
space. It is found that the maximum contribution to the variance comes from the islands 
immediately adjoining the central island. 



IV. POST NEWTONIAN WAVEFORM 

In the post-Newtonian case we have the five parameters, 

/X = /x" = /i\ /x^ /i^} = {A, 2nfsts, 27r/,ro, 27r/,ri}. 

Simulations (12000 realizations) were carried out again for an SNR of 10. The signal pa- 
rameters were tq = 5.558s and fi = 0.684s, corresponding to a binary comprised of a 
IOMq black hole and a IAMq neutron star. The simulation box taken was {5.058s< tq < 
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5.859s, 0.484s< ti < 0.985s}, with filter spacings of 10ms. in tq and 5ms. in ri. Our analysis 
of the post-Newtonian case will be very similar to the Newtonian case in section |T|. All the 
variables defined there can be defined for the post-Newtonian case also and we will use the 
same names for the variables to avoid using more symbols. 

Here again we diagonalise the covariance matrix by making a coordinate transformation 
from the n to the u coordinate system. The diagonal covariance matrix in the u parameters 
is given by. 



1 

l2 



0.0 0.0 0.0 0.0 
0.0 0.018 0.0 0.0 0.0 
0.0 0.0 1.1974 0.0 0.0 
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0.0 0.0 0.0 403.64 0.0 
yO.O 0.0 0.0 0.0 15601.0 
Therefore the root mean square errors in the parameters computed from the covariance 
matrix are, a^i = {.013,0.109,2.009, 12.49}. On the other hand the observed errors from 
the Monte Carlo simulations are E^, = {1.33574,7.43948,5.34089,23.049}. One finds as 
compared to the Newtonian case, the factors Hyi/a^i are larger on the average. The reason 
for this is that here we have one more parameter, which means there are more filters which 
can match with the data and consequently there are more outlying points. This is evident 
from Fig. |l^ as is explained below. 

In Figs. |l^ and we give the scatter plots of the four parameters i/', the former on the 



z/^ plane and the latter on the — z/^ plane. We observe that while there are gaps in 



distribution of the parameters and z/^ there are none in the distribution of the parameters 
z/^ and z/^. Here the phase parameter /i^ = $ is given by, 

$ = -0.109838z/^ - 0.621516z/^ + 0.695335z/^ + 0.343747z/^, 

and the variable Om is given by. 



Om = 4.1824lAz/^ + 0.892242Az/2 + 0.0186845Az/3 + 0.00272864Az/^ 
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Since both (Ji,z and cTyi are comparable, both these parameters will contribute to $ as 
opposed to the Newtonian case where only the parameter dominates. Thus we find that 
both the parameters and z/^ assume all their values in their respective ranges. The spacing 
between the islands in the z/^ and v"^ scatter plot is calculated to be 2tiGiICI} = 0.69 and 
27102/0^"^ = 3.90 respectively. This is consistent with Fig. |T^. In Fig. |T3| we notice that 



there are more islands to the right than there are to the left of the central island. This is 
caused by the simulation box which is asymmetrical, i.e. the parameters of the signal are 
not at the center of the simulation box. This is so since all combinations of tq and ti do not 
correspond to valid masses for the components of the binary. This leads to a shift in the 
mean of the estimated parameters from their actual values. The observed means in the Az/* 
are {0.10167, 0.574784, -0.282783, 1.56494}. For higher SNRs the estimated parameters will 
tend to lie only on the islands which lie close to the central island and consequently within 
a symmetrical region around the actual values of the parameters. Therefore the bias will 
disappear for the case of higher SNRs. 

In Fig. [T^l we plot the histogram of the distribution of the variable 6^ in the post 
Newtonian case. Here, for the same SNR of 10, there are more outlying points as compared 
to the Newtonian case. This is to be expected since we have an additional parameter, and 
there is an additional degree of freedom to find the filter which matches best with the signal. 
The variance now gets substantial contributions even from far away islands as opposed to 
the Newtonian case. 



In Fig. |T5| we plot the variables v/s 6^. Here again only r4 does not attain all its 
possible values within the central island, but now attains its entire range of values only 
when when we include two islands on either side. Moreover the overlap of values of 
in these islands is much more pronounced as opposed to the Newtonian case and further 
investigations are needed to identify the primary solution correctly. 

In Fig. we plot the probability distributions of the variables r,. The histograms 
are plotted using the Monte Carlo data and the continuous curve is a Gaussian with a 
variance of unity. It is clear that are Gaussian random variables to a good degree of 
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approximation. The approximation is better than in the Newtonian case (see Fig. P). As 
in the Newtonian case the rj and are related by only a constant factor. Therefore Ki are 
also Gaussian random variables. However, we note that this is true only for a special class 
of waveforms to which the chirp belongs. The chirp signal manifold is intrinsically flat and 
the parameterization of the waveform is such that the coordinates are essentially Cartesian. 
If some other parameterization is chosen {e.g. the chirp mass M), the coordinates will no 
longer be so and will not be Gaussian random variables even though the variables will 
remain Gaussian. 

We can again as in the Newtonian case, write down the expression for the pdf of the 
estimated parameters. We do not write it explicitly here since the expression is formally the 
same as in eqn. (p.29|) , except now the index a runs from to 4. We therefore only refer to 
eqn. ( p.29D for the required pdf. 



V. CONCLUSIONS 

In this paper we have addressed the question of wide discrepancies between the theoret- 
ically predicted lower bounds in the errors in the estimation of parameters and the errors 
observed in numerical experiments. We now summarize in this section, the main results of 
our paper and indicate future directions. 

• We find that the problem is of a global nature, in that the estimated values for the 
parameters fall into disjoint islands. Though there are very few points in the islands 
which are far from the actual value of the parameters, they contribute substantially 
to the variance. Thus the discrepancy between the Monte Carlo simulations and the 
Cramer Rao bounds cannot be resolved by using perturbative analysis. The restriction 
of the parameter $ to the range [— vr, vr] plays a major role in the non local distribution 
of parameters, as explained in the text. 

• The problem is more transparent when we reparameterize the chirp signal so that the 
covariance matrix is diagonal in the new parameters. The parameters v correspond to 
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choosing orthogonal coordinates on the chirp manifold. Since the covariances are zero 
for the variables the pdf is computationally simpler to handle. 

A statistical model has been given which matches well with the distribution obtained 
from Monte Carlo simulations. The model is derived from geometrical considerations. 
We have identified certain variables as Gaussian random variables and these play 
an important role in writing down the expression for the distribution of the estimated 
parameters. 

Since the distribution of the estimated parameters is multimodal, the variance is not a 
good indicator of the performance of the MLE. A more reasonable indicator would be 
to judge the performance of the MLE by means of confidence intervals i.e. compute 
the probability that the estimated parameters lie within a certain region around the 
actual value. As a concrete example, we compute the probability that the deviation of 
the estimated parameter is less than thrice the root mean square error at that SNR as 
predicted by the covariance matrix. We will use the symbol Pso- to denote the fraction 
of points which lie within the 3a range for a given parameter. However this criterion 
will in general be dependent on the specific parameter chosen. Since tq is one such 
physical parameter we use this to compute P^^. 



In Figures and |T8| we plot P30- v/s SNR for the parameter tq in the Newtonian and 
the post-Newtonian cases respectively. The values of P30- for the parameters such as 
To and Ti will be independent of the actual value of the chirp times. For this purpose 
we use the results of simulations carried out earlier in |jl9[. It is to be noted that 
whereas the simulations carried out in this paper use a single curve to fit the noise 
power spectrum, S'„(/), we had used a more accurate representation of S'„(/) in our 
earlier simulations |[T^ . We see that the MLE works moderately well even at low SNRs 
of p ~ 10. It is to be remarked that the assessment of an estimator depends upon how 
we use the estimates to calculate astrophysically interesting quantities. 
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• We required about 2 days of compution on a 300 MFlops (peak rating) machine to 
generate the results of this paper. The use of an integration routine specifically suited 
to the integrand, and/or the use of lookup tables for computing the integrand, would 
further speed up the computation substantially. The intrinsic flatness of the manifold 
turns out to be very convenient for our purpose. This property holds true for the 
first post-Newtonian waveform as well. There is one more parameter in the post- 
Newtonian waveform and consequently one more integration to perform to get the the 
marginal distribution of the parameters. For higher post-Newtonian corrections and/or 
for inclusion of parameters such as spins, it might be computationally expensive to 
compute the marginal distributions. However, it is to be noted that performing Monte 
Carlo simulations in such cases would also call for a huge amount of computational 
effort. A further research into the above issues is in progress. 
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FIGURES 

FIG. 1. Plot of /(/) v/s /. 

FIG. 2. Variation of root mean square errors in the parameters of the Newtonian waveform 
with SNR. 

FIG. 3. Plot of v/s /; Ai/^ = {1.62, -8.49, 7.3}. Om has a value close to -Air in the region 
around 135Hz, where /(/) attains its maximum. 

FIG. 4. Scatter plot of 9rn\ 12000 realizations of noise have been considered. 

FIG. 5. Scatter plot on the — plane 

FIG. 6. Scatter plot on the — z/^plane 

FIG. 7. Schematic illustration of the geometric picture discussed in the text. 

FIG. 8. Marginal Porbability distributions of the variables r*. 

FIG. 9. Plot of the variables v/s Qm- 

FIG. 10. Histogram of the variable d^- 

FIG. 11. Comparision between the statistical model and the Monte Carlo methods to arrive at 
the marginal distribution of the parameter v^. Part (c) and (d) zoom in to the first and second 
maxima to the right of the central maximum. 

FIG. 12. Histogram of 0^ for the post-Newtonian case. 
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FIG. 13. Scatter Plot on the - 



plane. 



FIG. 14. Scatter Plot on the - 



z/^ plane. 



FIG. 15. Plot of the variables v/s 



FIG. 16. Distributions of the variables r* for the post-Newtonian case are plotted. The his- 
togram is obtained from the Monte-Carlo simulations whereas the continuuous curve is a unit 
variance Gaussian curve. 

FIG. 17. Plot of Pso- for the parameter tq v/s SNR for the Newtonian case. 

FIG. 18. Plot of P'icr for the parameter tq v/s SNR for the post- Newtonian case. 
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